{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-12-09T15:53:09.502685Z",
     "start_time": "2022-12-09T15:53:08.653426Z"
    }
   },
   "outputs": [],
   "source": [
    "import open3d as o3d\n",
    "import open3d.core as o3c\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "import os\n",
    "import sys\n",
    "\n",
    "# Only needed for tutorial, monkey patches visualization\n",
    "sys.path.append(\"..\")\n",
    "import open3d_tutorial as o3dtut\n",
    "# Change to True if you want to interact with the visualization windows\n",
    "o3dtut.interactive = not \"CI\" in os.environ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PointCloud\n",
    "This tutorial demonstrates basic usage of a point cloud.\n",
    "\n",
    "## PointCloud creation\n",
    "The first part of the tutorial shows how to construct a point cloud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-12-09T15:53:09.507147Z",
     "start_time": "2022-12-09T15:53:09.503719Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create a empty point cloud on CPU.\n",
    "pcd = o3d.t.geometry.PointCloud()\n",
    "print(pcd, \"\\n\")\n",
    "\n",
    "# To create a point cloud on CUDA, specify the device.\n",
    "# pcd = o3d.t.geometry.PointCloud(o3c.Device(\"cuda:0\"))\n",
    "\n",
    "# Create a point cloud from open3d tensor with dtype of float32.\n",
    "pcd = o3d.t.geometry.PointCloud(o3c.Tensor([[0, 0, 0], [1, 1, 1]], o3c.float32))\n",
    "print(pcd, \"\\n\")\n",
    "\n",
    "# Create a point cloud from open3d tensor with dtype of float64.\n",
    "pcd = o3d.t.geometry.PointCloud(o3c.Tensor([[0, 0, 0], [1, 1, 1]], o3c.float64))\n",
    "print(pcd, \"\\n\")\n",
    "\n",
    "# Create a point cloud from numpy array. The array will be copied.\n",
    "pcd = o3d.t.geometry.PointCloud(\n",
    "    np.array([[0, 0, 0], [1, 1, 1]], dtype=np.float32))\n",
    "print(pcd, \"\\n\")\n",
    "\n",
    "# Create a point cloud from python list.\n",
    "pcd = o3d.t.geometry.PointCloud([[0., 0., 0.], [1., 1., 1.]])\n",
    "print(pcd, \"\\n\")\n",
    "\n",
    "# Error creation. The point cloud must have shape of (N, 3).\n",
    "try:\n",
    "    pcd = o3d.t.geometry.PointCloud(o3c.Tensor([0, 0, 0, 0], o3c.float32))\n",
    "except:\n",
    "    print(f\"Error creation. The point cloud must have shape of (N, 3).\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `PointCloud` can be created on both CPU and GPU, and with different data types. The device of the point cloud will be the same as the device of the input tensor.\n",
    "\n",
    "Besides, `PointCloud` can be also created by python dict with multiple attributes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-12-09T15:53:15.378492Z",
     "start_time": "2022-12-09T15:53:15.375715Z"
    }
   },
   "outputs": [],
   "source": [
    "map_to_tensors = {}\n",
    "\n",
    "# - The \"positions\" attribute must be specified.\n",
    "# - Common attributes include \"colors\" and \"normals\". \n",
    "# - You may also use custom attributes, such as \"labels\".\n",
    "# - The value of an attribute could be of any shape and dtype. Its correctness\n",
    "#   will only be checked when the attribute is used by some algorithms.\n",
    "map_to_tensors[\"positions\"] = o3c.Tensor([[0, 0, 0], [1, 1, 1]], o3c.float32)\n",
    "map_to_tensors[\"normals\"] = o3c.Tensor([[0, 0, 1], [0, 0, 1]], o3c.float32)\n",
    "map_to_tensors[\"labels\"] = o3c.Tensor([0, 1], o3c.int64)\n",
    "\n",
    "pcd = o3d.t.geometry.PointCloud(map_to_tensors)\n",
    "print(pcd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point cloud attributes setter and getter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-12-09T15:53:24.681751Z",
     "start_time": "2022-12-09T15:53:24.678427Z"
    }
   },
   "outputs": [],
   "source": [
    "pcd = o3d.t.geometry.PointCloud(o3c.Tensor([[0, 0, 0], [1, 1, 1]], o3c.float32))\n",
    "# Set attributes.\n",
    "pcd.point.normals = o3c.Tensor([[0, 0, 1], [0, 0, 1]], o3c.float32)\n",
    "pcd.point.colors = o3c.Tensor([[1, 0, 0], [0, 1, 0]], o3c.float32)\n",
    "pcd.point.labels = o3c.Tensor([0, 1], o3c.int64)\n",
    "print(pcd, \"\\n\")\n",
    "\n",
    "# Set by numpy array or python list.\n",
    "pcd.point.normals = np.array([[0, 0, 1], [0, 0, 1]], dtype=np.float32)\n",
    "pcd.point.intensity = [0.4, 0.4]\n",
    "print(pcd, \"\\n\")\n",
    "\n",
    "# Get attributes.\n",
    "posisions = pcd.point.positions\n",
    "print(\"posisions: \")\n",
    "print(posisions, \"\\n\")\n",
    "labels = pcd.point.labels\n",
    "print(\"labels: \")\n",
    "print(labels, \"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conversion between tensor and legacy point cloud\n",
    "PointCloud can be converted to/from legacy `open3d.geometry.PointCloud`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "legacy_pcd = pcd.to_legacy()\n",
    "print(legacy_pcd, \"\\n\")\n",
    "\n",
    "tensor_pcd = o3d.t.geometry.PointCloud.from_legacy(legacy_pcd)\n",
    "print(tensor_pcd, \"\\n\")\n",
    "\n",
    "# Convert from legacy point cloud with data type of float64.\n",
    "tensor_pcd_f64 = o3d.t.geometry.PointCloud.from_legacy(legacy_pcd, o3c.float64)\n",
    "print(tensor_pcd_f64, \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize point cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Load a ply point cloud, print it, and render it\")\n",
    "ply_point_cloud = o3d.data.PLYPointCloud()\n",
    "pcd = o3d.t.io.read_point_cloud(ply_point_cloud.path)\n",
    "print(pcd)\n",
    "o3d.visualization.draw_geometries([pcd.to_legacy()],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`read_point_cloud` reads a point cloud from a file. It tries to decode the file based on the extension name. For a list of supported file types, refer to [File IO](file_io.ipynb).\n",
    "\n",
    "`draw_geometries` visualizes the point cloud. Use a mouse/trackpad to see the geometry from different view points.\n",
    "\n",
    "It looks like a dense surface, but it is actually a point cloud rendered as surfels. The GUI supports various keyboard functions. For instance, the `-` key reduces the size of the points (surfels).\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "Press the `H` key to print out a complete list of keyboard instructions for the GUI. For more information of the visualization GUI, refer to [Visualization](visualization.ipynb) and [Customized visualization](../visualization/customized_visualization.rst).\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "On macOS, the GUI window may not receive keyboard events. In this case, try to launch Python with `pythonw` instead of `python`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downsampling\n",
    "This section provides several downsamle methods of a point cloud.\n",
    "### Voxel downsampling\n",
    "Voxel downsampling uses a regular voxel grid to create a uniformly downsampled point cloud from an input point cloud. It is often used as a pre-processing step for many point cloud processing tasks. The algorithm operates in two steps:\n",
    "\n",
    "1. Points are bucketed into voxels.\n",
    "2. Each occupied voxel generates exactly one point by averaging all points inside. \n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "Currently, the method returns the voxel coordinates of the point cloud.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Downsample the point cloud with a voxel of 0.03\")\n",
    "downpcd = pcd.voxel_down_sample(voxel_size=0.03)\n",
    "o3d.visualization.draw_geometries([downpcd.to_legacy()],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Farthest point downsampling\n",
    "Farthest point sampling samples the point cloud by selecting the farthest point from the current selected point iteratively. It is used to sample the point cloud to a fixed number of points which holds the maximum geometrical information of the original point cloud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Downsample the point cloud by selecting 5000 farthest points.\")\n",
    "downpcd_farthest = pcd.farthest_point_down_sample(5000)\n",
    "o3d.visualization.draw_geometries([downpcd_farthest.to_legacy()],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open3D also provides `uniform_down_sample` and `random_down_sample` for point cloud downsampling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vertex normal estimation\n",
    "Another basic operation for point cloud is normal estimation.\n",
    "Press `N` to see point normals. The keys `-` and `+` can be used to control the length of the normal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Recompute the normal of the downsampled point cloud using hybrid nearest neighbor search with 30 max_nn and radius of 0.1m.\")\n",
    "downpcd.estimate_normals(max_nn=30, radius=0.1)\n",
    "o3d.visualization.draw_geometries([downpcd.to_legacy()],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024],\n",
    "                                  point_show_normal=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`estimate_normals` computes the normal for every point. The function finds adjacent points and calculates the principal axis of the adjacent points using covariance analysis.\n",
    "\n",
    "The two key arguments `radius = 0.1` and `max_nn = 30` specifies search radius and maximum nearest neighbor. It has 10cm of search radius, and only considers up to 30 neighbors to save computation time. If only `max_nn` or `radius` is specified (`downpcd.estimate_normals(30, None)` or `downpcd.estimate_normals(None, 0.01)`), the function will use knn or radius search respectively. \n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "It is always recommended to specify both `radius` and `max_nn`, especially if the point cloud is on CUDA device.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "The covariance analysis algorithm produces two opposite directions as normal candidates. Without knowing the global structure of the geometry, both can be correct. This is known as the normal orientation problem. Open3D tries to orient the normal to align with the original normal if it exists. Otherwise, Open3D does a random guess. Further orientation functions such as `orient_normals_to_align_with_direction` and `orient_normals_towards_camera_location` need to be called if the orientation is a concern.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Access estimated normals\n",
    "The estimated normals can be access and converted to numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normals = downpcd.point.normals\n",
    "print(\"Print first 5 normals of the downsampled point cloud.\")\n",
    "print(normals[:5], \"\\n\")\n",
    "print(\"Convert normals tensor into numpy array.\")\n",
    "normals_np = normals.numpy()\n",
    "print(normals_np[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Crop point cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Load a polygon volume and use it to crop the original point cloud\")\n",
    "demo_crop_data = o3d.data.DemoCropPointCloud()\n",
    "pcd = o3d.t.io.read_point_cloud(demo_crop_data.point_cloud_path)\n",
    "vol = o3d.visualization.read_selection_polygon_volume(demo_crop_data.cropped_json_path)\n",
    "chair = vol.crop_point_cloud(pcd.to_legacy())\n",
    "o3d.visualization.draw_geometries([chair],\n",
    "                                  zoom=0.7,\n",
    "                                  front=[0.5439, -0.2333, -0.8060],\n",
    "                                  lookat=[2.4615, 2.1331, 1.338],\n",
    "                                  up=[-0.1781, -0.9708, 0.1608])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`read_selection_polygon_volume` reads a json file that specifies polygon selection area. `vol.crop_point_cloud(pcd)` filters out points. Only the chair remains.\n",
    "\n",
    "We can also get the cropped indices of the point cloud using `crop_in_polygon`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indices = vol.crop_in_polygon(pcd.to_legacy())\n",
    "print(f\"Cropped indices length: {len(indices)}\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paint point cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Paint point cloud.\")\n",
    "pcd.paint_uniform_color([1, 0.706, 0])\n",
    "o3d.visualization.draw_geometries([pcd.to_legacy()],\n",
    "                                  zoom=0.7,\n",
    "                                  front=[0.5439, -0.2333, -0.8060],\n",
    "                                  lookat=[2.4615, 2.1331, 1.338],\n",
    "                                  up=[-0.1781, -0.9708, 0.1608])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`paint_uniform_color` paints all the points to a uniform color. The color is in RGB space, [0, 1] range."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bounding volumes\n",
    "The `PointCloud` geometry type has bounding volumes as all other geometry types in Open3D. Currently, Open3D implements an `AxisAlignedBoundingBox` and an `OrientedBoundingBox` that can also be used to crop the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aabb = pcd.get_axis_aligned_bounding_box()\n",
    "aabb.set_color(o3c.Tensor([1, 0, 0], o3c.float32))\n",
    "obb = pcd.get_oriented_bounding_box()\n",
    "obb.set_color(o3c.Tensor([0, 1, 0], o3c.float32))\n",
    "o3d.visualization.draw_geometries(\n",
    "    [pcd.to_legacy(), aabb.to_legacy(),\n",
    "     obb.to_legacy()],\n",
    "    zoom=0.7,\n",
    "    front=[0.5439, -0.2333, -0.8060],\n",
    "    lookat=[2.4615, 2.1331, 1.338],\n",
    "    up=[-0.1781, -0.9708, 0.1608])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point cloud outlier removal\n",
    "When collecting data from scanning devices, the resulting point cloud tends to contain noise and artifacts that one would like to remove. This demo below addresses the outlier removal features of Open3D.\n",
    "\n",
    "### Prepare input data\n",
    "A point cloud is loaded and downsampled using `voxel_downsample`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Load a ply point cloud, print it, and render it\")\n",
    "sample_pcd_data = o3d.data.PCDPointCloud()\n",
    "pcd = o3d.t.io.read_point_cloud(sample_pcd_data.path)\n",
    "o3d.visualization.draw_geometries([pcd.to_legacy()],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])\n",
    "\n",
    "print(\"Downsample the point cloud with a voxel of 0.02\")\n",
    "voxel_down_pcd = pcd.voxel_down_sample(voxel_size=0.02)\n",
    "o3d.visualization.draw_geometries([voxel_down_pcd.to_legacy()],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, use `uniform_down_sample` to downsample the point cloud by collecting every n-th points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Every 5th points are selected\")\n",
    "uni_down_pcd = pcd.uniform_down_sample(every_k_points=5)\n",
    "o3d.visualization.draw_geometries([uni_down_pcd.to_legacy()],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select down sample\n",
    "The following helper function uses `select_by_mask`, which takes a binary mask to output only the selected points. The selected points and the non-selected points are visualized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def display_inlier_outlier(cloud : o3d.t.geometry.PointCloud, mask : o3c.Tensor):\n",
    "    inlier_cloud = cloud.select_by_mask(mask)\n",
    "    outlier_cloud = cloud.select_by_mask(mask, invert=True)\n",
    "\n",
    "    print(\"Showing outliers (red) and inliers (gray): \")\n",
    "    outlier_cloud = outlier_cloud.paint_uniform_color([1.0, 0, 0])\n",
    "    inlier_cloud.paint_uniform_color([0.8, 0.8, 0.8])\n",
    "    inlier_cloud = o3d.visualization.draw_geometries([inlier_cloud.to_legacy(), outlier_cloud.to_legacy()],\n",
    "                                      zoom=0.3412,\n",
    "                                      front=[0.4257, -0.2125, -0.8795],\n",
    "                                      lookat=[2.6172, 2.0475, 1.532],\n",
    "                                      up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical outlier removal\n",
    "`statistical_outlier_removal` removes points that are further away from their neighbors compared to the average for the point cloud. It takes two input parameters:\n",
    "\n",
    "- `nb_neighbors`, which specifies how many neighbors are taken into account in order to calculate the average distance for a given point.\n",
    "- `std_ratio`, which allows setting the threshold level based on the standard deviation of the average distances across the point cloud. The lower this number the more aggressive the filter will be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Statistical oulier removal\")\n",
    "cl, ind = voxel_down_pcd.remove_statistical_outliers(nb_neighbors=20,\n",
    "                                                     std_ratio=2.0)\n",
    "display_inlier_outlier(voxel_down_pcd, ind)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Radius outlier removal\n",
    "`radius_outlier_removal` removes points that have few neighbors in a given sphere around them. Two parameters can be used to tune the filter to your data:\n",
    "\n",
    "- `nb_points`, which lets you pick the minimum amount of points that the sphere should contain.\n",
    "- `radius`, which defines the radius of the sphere that will be used for counting the neighbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Radius oulier removal\")\n",
    "cl, ind = voxel_down_pcd.remove_radius_outliers(nb_points=16, search_radius=0.05)\n",
    "display_inlier_outlier(voxel_down_pcd, ind)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convex hull\n",
    "The convex hull of a point cloud is the smallest convex set that contains all points. Open3D contains the method `compute_convex_hull` that computes the convex hull of a point cloud. The implementation is based on [Qhull](http://www.qhull.org/).\n",
    "\n",
    "In the example code below we compute the convex hull that is returned as a triangle mesh. Then, we visualize the convex hull as a red `LineSet`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bunny = o3d.data.BunnyMesh()\n",
    "pcd = o3d.t.io.read_point_cloud(bunny.path)\n",
    "\n",
    "hull = pcd.compute_convex_hull()\n",
    "hull_ls = o3d.geometry.LineSet.create_from_triangle_mesh(hull.to_legacy())\n",
    "hull_ls.paint_uniform_color((1, 0, 0))\n",
    "o3d.visualization.draw_geometries([pcd.to_legacy(), hull_ls])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DBSCAN clustering\n",
    "Given a point cloud from e.g. a depth sensor we want to group local point cloud clusters together. For this purpose, we can use clustering algorithms. Open3D implements DBSCAN [\\[Ester1996\\]](../reference.html#Ester1996) that is a density based clustering algorithm. The algorithm is implemented in `cluster_dbscan` and requires two parameters: `eps` defines the distance to neighbors in a cluster and `min_points` defines the minimum number of points required to form a cluster. The function returns `labels`, where the label `-1` indicates noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ply_point_cloud = o3d.data.PLYPointCloud()\n",
    "pcd = o3d.t.io.read_point_cloud(ply_point_cloud.path)\n",
    "\n",
    "with o3d.utility.VerbosityContextManager(\n",
    "        o3d.utility.VerbosityLevel.Debug) as cm:\n",
    "    labels = pcd.cluster_dbscan(eps=0.02, min_points=10, print_progress=True)\n",
    "\n",
    "max_label = labels.max().item()\n",
    "print(f\"point cloud has {max_label + 1} clusters\")\n",
    "colors = plt.get_cmap(\"tab20\")(\n",
    "        labels.numpy() / (max_label if max_label > 0 else 1))\n",
    "colors = o3c.Tensor(colors[:, :3], o3c.float32)\n",
    "colors[labels < 0] = 0\n",
    "pcd.point.colors = colors\n",
    "o3d.visualization.draw_geometries([pcd.to_legacy()],\n",
    "                                  zoom=0.455,\n",
    "                                  front=[-0.4999, -0.1659, -0.8499],\n",
    "                                  lookat=[2.1813, 2.0619, 2.0999],\n",
    "                                  up=[0.1204, -0.9852, 0.1215])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "This algorithm precomputes all neighbors in the epsilon radius for all points. This can require a lot of memory if the chosen epsilon is too large.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plane segmentation\n",
    "Open3D also supports segmententation of geometric primitives from point clouds using RANSAC. To find the plane with the largest support in the point cloud, we can use `segment_plane`. The method has four arguments: `distance_threshold` defines the maximum distance a point can have to an estimated plane to be considered an inlier, `ransac_n` defines the number of points that are randomly sampled to estimate a plane, `num_iterations` defines how often a random plane is sampled and verified, and `probability` defined the expected probability of finding the optimal plane. The function then returns the plane as $(a,b,c,d)$ such that for each point $(x,y,z)$ on the plane we have $ax + by + cz + d = 0$. The function further returns a indices of the inlier points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_pcd_data = o3d.data.PCDPointCloud()\n",
    "pcd = o3d.t.io.read_point_cloud(sample_pcd_data.path)\n",
    "plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,\n",
    "                                            ransac_n=3,\n",
    "                                            num_iterations=1000, \n",
    "                                            probability=0.9999)\n",
    "[a, b, c, d] = plane_model.numpy().tolist()\n",
    "print(f\"Plane equation: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0\")\n",
    "\n",
    "inlier_cloud = pcd.select_by_index(inliers)\n",
    "inlier_cloud = inlier_cloud.paint_uniform_color([1.0, 0, 0])\n",
    "outlier_cloud = pcd.select_by_index(inliers, invert=True)\n",
    "o3d.visualization.draw_geometries([inlier_cloud.to_legacy(), outlier_cloud.to_legacy()],\n",
    "                                  zoom=0.8,\n",
    "                                  front=[-0.4999, -0.1659, -0.8499],\n",
    "                                  lookat=[2.1813, 2.0619, 2.0999],\n",
    "                                  up=[0.1204, -0.9852, 0.1215])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To achieve the stable results with specific random seed, you can set `probability` to 1.0, which will force the executed iteration equal to `num_iterations`. (Since `0.16.0`, the `segment_plane` function is parallel, and the iteration will be updated by the equation: `iter = log(1 - probability) / log(1 - fitness ^ ransac_n)`, where `fitness` is the ratio of inliers number and total number of points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "o3d.utility.random.seed(0)\n",
    "\n",
    "plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,\n",
    "                                            ransac_n=3,\n",
    "                                            num_iterations=1000, \n",
    "                                            probability=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hidden point removal\n",
    "Imagine you want to render a point cloud from a given view point, but points from the background leak into the foreground because they are not occluded by other points. For this purpose we can apply a hidden point removal algorithm. In Open3D the method by [\\[Katz2007\\]](../reference.html#Katz2007) is implemented that approximates the visibility of a point cloud from a given view without surface reconstruction or normal estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Convert mesh to a point cloud and estimate dimensions\")\n",
    "armadillo = o3d.data.ArmadilloMesh()\n",
    "mesh = o3d.io.read_triangle_mesh(armadillo.path)\n",
    "# Tensor TriangleMesh not supported this function yet.\n",
    "mesh.compute_vertex_normals()\n",
    "\n",
    "pcd = mesh.sample_points_poisson_disk(5000)\n",
    "diameter = np.linalg.norm(\n",
    "    np.asarray(pcd.get_max_bound()) - np.asarray(pcd.get_min_bound()))\n",
    "o3d.visualization.draw_geometries([pcd])\n",
    "\n",
    "print(\"Define parameters used for hidden_point_removal\")\n",
    "camera = o3d.core.Tensor([0, 0, diameter], o3d.core.float32)\n",
    "radius = diameter * 100\n",
    "\n",
    "print(\"Get all points that are visible from given view point\")\n",
    "pcd = o3d.t.geometry.PointCloud.from_legacy(pcd)\n",
    "_, pt_map = pcd.hidden_point_removal(camera, radius)\n",
    "pcd = pcd.select_by_index(pt_map)\n",
    "\n",
    "print(\"Visualize result\")\n",
    "o3d.visualization.draw_geometries([pcd.to_legacy()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boundary detection\n",
    "Open3D implements the boundary detection algorithm inspired by [PCL](https://pointclouds.org/documentation/classpcl_1_1_boundary_estimation.html). The algorithm find the boundary points among a unordered point cloud by analyzing the angle among the normals of a point and its neighbors. The method has three arguments: `radius` and `max_nn` specify the hybrid nearest search parameters; The `angle_threshold` defines the maximum angle between the normals of a point and its neighbors to be considered as a boundary point. The functions returns the boundary points and a boolean tensor of the same size as the input point cloud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ply_point_cloud = o3d.data.DemoCropPointCloud()\n",
    "pcd = o3d.t.io.read_point_cloud(ply_point_cloud.point_cloud_path)\n",
    "\n",
    "boundarys, mask = pcd.compute_boundary_points(0.02, 30)\n",
    "# TODO: not good to get size of points.\n",
    "print(f\"Detect {boundarys.point.positions.shape[0]} bnoundary points from {pcd.point.positions.shape[0]} points.\")\n",
    "\n",
    "boundarys = boundarys.paint_uniform_color([1.0, 0.0, 0.0])\n",
    "pcd = pcd.paint_uniform_color([0.6, 0.6, 0.6])\n",
    "o3d.visualization.draw_geometries([pcd.to_legacy(), boundarys.to_legacy()],\n",
    "                                      zoom=0.3412,\n",
    "                                      front=[0.3257, -0.2125, -0.8795],\n",
    "                                      lookat=[2.6172, 2.0475, 1.532],\n",
    "                                      up=[-0.0694, -0.9768, 0.2024])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "c36d4cd2cfce5c26977a9560a8c96efbe52e59f8182dcd6a297957684d25bbef"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
